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1.  Introduction 

The  numerical  solution  of  problems  in  very  large  domains  has  been  an  active  area  of  research  for  the  last  three  decades 
[1].  In  the  context  of  fields  of  application  like  acoustics,  electromagnetics,  meteorology,  solid  geophysics  and  aerodynamics 
where  dispersive  wave  problems  arise,  there  are  several  available  methods.  The  use  of  non-reflecting  boundary  conditions 
(NRBCs)  is  one  such  method. 

The  method  of  NRBCs  can  be  described  as  follows.  First,  the  infinite  domain  is  truncated  via  an  artificial  boundary  B,  thus 
dividing  the  original  domain  into  a  finite  computational  domain  Q  and  a  residual  infinite  domain  D.  Then,  a  special  boundary 
condition  is  imposed  on  B  in  order  to  complete  the  statement  of  the  problem  in  Q  (i.e.,  make  the  solution  in  Q  unique)  and, 
most  importantly,  to  ensure  that  little  or  no  spurious  wave  reflections  occur  from  B.  This  boundary  condition  is  called  a 
NRBC,  although  other  names  are  often  used  (see  [2]).  Finally,  the  problem  is  solved  numerically  in  Q ,  say  by  the  spectral 
element  (SE)  method.  The  SE  method  is  a  generalized  high-order  finite  element  method  where  the  integration  and  interpo¬ 
lation  points  are  selected  carefully  in  order  to  yield  accurate  but  efficient  solutions.  See  Patera  [3]  and  Giraldo-Restelli  [4]  for 
more  details  on  this  method. 

Fig.  1  illustrates  the  NRBC  set-up  using  an  infinite  wave-guide.  Here,  the  artificial  boundary  B  extends  from  the  southern 
(Is)  to  the  northern  (rN)  boundaries  of  the  wave-guide,  thus  creating  the  east  (rE)  and  west  (rw)  boundaries  of  Q  at 
x  =  xE,xw  respectively.  Outside  of  the  area  enclosed  by  these  boundaries  is  the  residual  infinite  domain  D. 

Naturally,  the  quality  of  the  numerical  solution  strongly  depends  on  the  properties  of  the  NRBC  employed.  In  fact,  it  can 
be  argued  that  for  many  applications  (e.g.,  limited-area  atmospheric  modeling  [4])  the  quality  of  the  solution  is  almost 
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Fig.  1.  An  infinite  channel  domain  Q  truncated  by  artificial  boundaries  rw  and  rE. 


entirely  due  to  the  accuracy  of  the  NRBC.  In  the  last  35  years  or  so,  much  research  has  been  done  to  develop  NRBCs  that,  after 
discretization,  lead  to  a  scheme  that  is  stable,  accurate,  efficient  and  easy  to  implement.  See  [5-8]  for  recent  reviews  on  the 
subject.  Of  course,  it  is  difficult  to  find  a  single  NRBC  which  is  ideal  in  all  respects  and  all  cases;  this  is  why  the  quest  for 
better  NRBCs  and  their  associated  discretization  schemes  continues. 

Recently,  high-order  local  NRBCs  have  been  introduced.  Sequences  of  increasing-order  NRBCs  have  been  available  for  a 
long  time  (e.g.,  the  Bayliss-Turkel  conditions  [9]  constitute  such  a  sequence),  but  they  had  been  regarded  as  impractical  be¬ 
yond  2nd  or  3rd  order  from  the  implementation  point  of  view.  Only  since  the  mid  1990s  have  practical  high-order  NRBCs 
been  devised.  The  present  paper  is  concerned  with  such  a  high-order  NRBC  scheme. 

The  first  high-order  local  NRBC  was  proposed  by  Collino  [10],  for  two-dimensional  time-dependent  waves  in  rectan¬ 
gular  domains.  Its  construction  requires  the  solution  of  the  one-dimensional  wave  equation  on  8.  Grote  and  Keller  [11] 
developed  a  high-order  converging  NRBC  for  the  three-dimensional  time-dependent  wave  equation,  based  on  spherical 
harmonic  transformations.  Sofronov  [12,13]  proposed  exact  boundary  conditions  for  the  three-  and  two-dimensional 
wave  equations  in  spherical  and  polar  coordinates,  respectively  (it  is  proved  that  NRBCs  demonstrated  in  [11,12]  are  re¬ 
duced  to  each  other).  Hagstrom  and  Hariharan  [14]  constructed  high-order  NRBCs  for  the  two-  and  three-dimensional 
time-dependent  wave  equations  based  on  the  analytic  series  representation  for  the  outgoing  solutions  of  these  equa¬ 
tions.  For  time-dependent  waves  in  a  two-dimensional  wave-guide,  Guddati  and  Tassoulas  [15]  devised  a  high-order 
NRBC  by  using  rational  approximations  and  recursive  continued  fractions.  Givoli  [16]  has  shown  how  to  derive  high-or¬ 
der  NRBCs  for  a  general  class  of  wave  problems,  leading  to  a  symmetric  FE  formulation.  In  [17],  this  methodology  was 
applied  to  the  particular  case  of  time-harmonic  waves,  using  optimally  localized  Dirichlet-to-Neumann  (DtN)  NRBCs  (see 
also  [18]). 

The  starting  point  for  the  family  of  NRBCs  discussed  here  is  the  condition  devised  by  Higdon  [19],  which  was  designed  for 
low-order  finite  difference  schemes.  Givoli  and  Neta  [20]  directly  extended  the  Higdon  scheme  to  high-order  finite  difference 
discretizations.  They  later  extended  this  formulation  to  one  that  does  not  involve  any  high  derivatives  (hereafter  referred  to 
as  the  G-N  formulation).  The  elimination  of  all  high-order  derivatives  is  enabled  through  the  introduction  of  special  auxiliary 
variables  on  8.  This  construction  demonstrated  in  [21,22]  for  finite  differences  was  further  extended  in  [23]  for  finite  ele¬ 
ment  schemes  to  solve  the  dispersive  wave  equation.  Hagstrom  and  Warburton  [24]  used  the  Higdon  and  auxiliary  variable 
framework  to  develop  a  symmetric  boundary  formulation  in  a  full-space  configuration  where  special  corner  compatibility 
conditions  were  developed  for  the  non-dispersive  wave  equation.  Extensions  and  comparisons  between  the  two  methods 
were  published  by  Givoli  and  coworkers  in  [25,26]. 

In  the  present  paper,  the  G-N  auxiliary  formulation  of  the  Higdon  NRBC  is  implemented  on  a  reduced  form  of  the  line¬ 
arized  shallow  water  equations.  This  formulation  has  been  previously  demonstrated  in  both  finite  difference  and  finite  ele¬ 
ment  schemes  to  arbitrarily  high  NRBC  order,  however,  accuracy  gains  realized  by  increasing  the  NRBC  order  slowed 
significantly  after  orders  2  or  3.  The  formulation  used  here  seeks  to  remedy  this  limitation  by  using  a  high-order  treatment 
of  space  (SE)  and  time  (Runge-Kutta)  to  show  the  benefits  of  using  a  high-order  boundary  (G-N)  scheme.  Specifically,  the 
interior  and  boundary  formulations  are  discretized  using  high-order  basis  functions  in  a  stable,  equal-order  interpolation 
scheme  for  all  the  variables  (this  does  not  violate  the  inf-sup  condition).  High-order  time  integration  is  performed  as  well 
in  an  effort  to  balance  all  of  the  errors  involved  with  the  numerical  solution.  The  computational  effort  associated  with  the 
high-order  boundary  scheme  can  be  shown  to  grow  only  linearly  with  the  order. 

It  should  be  noted  that  the  only  other  spectral  element,  high-order  boundary  approach  (to  the  authors’  knowledge)  is  [27] 
that  is  in  press.  That  paper  shows  similar  results  for  the  non-dispersive  wave  equation  on  a  semi-infinite  channel  using  the 
Hagstrom-Warburton  formulation.  The  key  difference  in  our  work  is  that  we  use  high-order  space,  boundary  and  time  inte¬ 
gration  in  both  a  dispersive  and  non-dispersive  wave  equation  setting. 

The  following  is  the  outline  of  the  rest  of  this  paper:  In  Section  2,  the  problem  under  investigation  is  stated.  In  Sec¬ 
tion  3,  an  overview  of  the  G-N  auxiliary  formulation  is  presented.  In  Section  4,  a  SE  semi-discrete  formulation  is  con¬ 
structed  which  incorporates  a  NRBC  of  this  family  with  any  desired  order.  In  Section  5  the  time-integrator  used  to 
march  the  equations  in  time  is  discussed.  The  performance  of  the  method  is  demonstrated  in  Section  6  via  numerical 
examples. 
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2.  Statement  of  the  problem 

As  a  model  serving  for  introducing  the  ideas  developed  here,  the  shallow  water  equations  are  considered.  To  fix  ideas, 
some  specific  equations  and  boundary  conditions  are  chosen  here.  In  the  channel  (Fig.  1)  the  2-D  shallow  water  equations 
(SWEs)  are 

dtu  +  liudxu  +  livdyU  -fv  =  - gdxrj , 

dtv  +  imdx  v  +  /ivdyV+fu  =  -gdyf],  (1) 

dtrj  +  iiudxr]  +  jivdyf]  +  (h0  +  m)(<9*u  +  dyv)  =  0. 

Here  t  is  time,  u(x,y,  t)  and  v(x,y,  t)  are  the  unknown  velocities  in  the  x  and  y  directions,  h0  is  the  given  water  layer  thickness 
(in  the  direction  normal  to  the  xy  plane),  rj(x,y,  t)  is  the  unknown  water  elevation  above  h0,/ is  the  Coriolis  parameter,  and  g 
is  the  gravity  acceleration.  We  use  the  following  shorthand  for  partial  derivatives: 


The  parameter  ^  is  1  for  the  nonlinear  SWEs,  and  is  0  for  the  linearized  SWEs  with  vanishing  mean  flow.  We  shall  consider 
the  latter  as  a  special  case  in  this  analysis.  It  can  be  shown  that  this  system  of  equations  can  be  reduced  to  the  dispersive 
wave  equation,  also  called  Klein-Gordon  equation  (KGE): 

d2tr,-C20V2r,+f2r,  =  S.  (2) 

where  C2Q  =  gh0.  The  goal  of  this  analysis  is  to  solve  the  analogous  homogeneous  problem 

d2u  -  C20V2u  +f2u  =  0.  (3) 

in  the  finite  domain  Q  via  SEs.  On  the  south  and  north  channel  walls  ( rs  and  rN)  we  have  no  normal  flow,  i.e: 

|]  =  0  on  rN  and  rs.  (4) 

At  x  — >  ±oo  the  solution  is  known  to  be  bounded  and  not  to  include  any  incoming  waves.  To  complete  the  statement  of  the 
problem,  initial  conditions  must  be  given  in  the  entire  domain  Q  u  D.  To  this  end,  the  artificial  boundary  B  is  now  introduced 
at  x  =  xE  and  x  =  xw  to  divide  between  Q  and  D  as  shown  in  Fig.  1. 


3.  Auxiliary  variables 


We  present  a  brief  summary  of  the  G-N  auxiliary  variable  process  as  described  in  [23].  This  auxiliary  formulation  begins 
with  the  Higdon  [19]  boundary  condition  given  by: 


Hj: 


u  =  0  on  rE. 


(5) 


Next,  we  introduce  the  auxiliary  functions  fa , . . . ,  which  are  defined  on  rE  and  rw  as  well  as  in  the  exterior  domain  D 
(see  Fig.  1).  Eventually  we  shall  use  these  functions  only  on  rE  and  rw,  but  the  derivation  requires  that  they  be  defined  in  D 
as  well,  or  at  least  in  a  non-vanishing  region  adjacent  to  rE  and  rw.  The  functions  fa  are  defined  via  the  relations 


(&+Tat) 

(5x+^5t) 


u  =  h, 
01  =  02 1 


(6) 


,=0. 

By  definition,  these  relations  hold  in  D,  and  also  on  rE  and  rw.  It  is  easy  to  see  that  (6),  when  imposed  as  boundary  condi¬ 
tions  on  rE  and  rw,  are  equivalent  to  the  single  boundary  condition  (5).  If  we  also  define 

0o  =  ui  0;  =  0,  (7) 


then  we  can  write  (6)  concisely  as 
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This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appearance  of  the  x-derivative  in  (8),  one 
cannot  discretize  the  on  the  boundary  rE  alone.  Therefore  we  shall  manipulate  (8)  in  order  to  get  rid  of  the  x-derivative. 

The  function  u  satisfies  the  dispersive  wave  Eq.  (3)  in  D.  It  is  easy  to  show  that  the  functions  satisfy  an  equation  like  (3), 
namely, 


d2x(f>j  +  ~  -  4 =  0- 


Here  we  need  the  assumption  that  C0  and  /  do  not  depend  on  x  or  on  t. 
Givoli  et  al.  [23]  have  shown  that  (8),  for  j  =  1, ...  ,J  -  1,  is  equivalent  to 


(9) 


1  1 


92 
dy 2 


(10) 


As  desired,  the  new  boundary  condition  (10)  does  not  involve  x-derivatives.  In  addition,  there  are  no  high  y-  and  t-derivatives 
in  (10)  beyond  second  order. 

Rewriting  (6),  (10)  and  (7),  the  new  formulation  of  the  Jth-order  NRBC  on  B  can  be  summarized  as  follows: 


„  ■  du 

Pou  +  Qj-h, 

(11) 

Pj'Pj  —  zJh-t  —  4>j- 1  + 

i  —  4>j+ 1  J  —  ft 

•  J- 1, 

(12) 

1  1  .  1 

J  C2  C2>  ^0  Ci- 

0  1  1 

ft  ~  c  +  C  ’ 

LJ  LJ+ 1 

II 

(13) 

d 

III 

3 

III 

o 

-e- 

(14) 

In  (12)  and  elsewhere,  a  dot  indicates  time  derivative  and  a  prime  indicates  differentiation  with  respect  toy  along  B,  i.e.,  the 
tangential  derivative  on  B. 


4.  Spectral  element  method 

As  a  starting  point,  the  semi-infinite  channel  is  considered.  This  set-up  is  identical  to  that  described  previously,  except 
that  at  rWt  a  Dirichlet  boundary  condition  u(0,y,  t)  =  0  is  prescribed;  see  Fig.  2. 

Now,  the  weak  form  of  the  problem  in  Q  is  constructed.  The  solution  is  sought  in  the  space  of  trial  functions, 

V  =  {u|u  g  H'(Q)  and  u  =  0  on  rw}.  (15) 

Here,  H1  (Q)  is  the  Sobolev  space  of  square-integrable  functions  and  their  first  derivatives  in  Q.  Now,  Eq.  (3)  is  multiplied  by 
the  basis  functions  xPi(x,y)  e  V  and  integrated  over  Q  so  the  weak  form  is  (after  integration  by  parts) 

[  WiUtt dQ  +  Cg  [  V'FiVudQ-Cl  [  'Ffi  Vudr+f2  [  %udQ  =  0.  (16) 

J  Q  J  Q  Jr  J Q 

Along  the  eastern  boundary,  the  weak  form  is  found  similarly  by  multiplying  (12)  by  the  test  function  t j  and  integrating  it 
over  rE.  After  integration  by  parts,  simplifying  and  using  (4)  to  deduce  that  the  boundary  term  vanishes,  this  yields: 

ft  [  dr E  —  (%j  I  drE  +  f  'Zjfij-]  dr E  +  a  I  dr E  =  f  ij^-,  dr e  (17) 

JrE  JrE  JrE  JrE  JrE 

for  j  =  1 , . . .  J  -  1 ,  ^  6  H1  (Te)  and  any  Tj  e  H1  (T£). 


VN 


rv 

b  \  rE 

u 

D 

Fig.  2.  A  semi-infinite  channel  domain  Q  truncated  by  artificial  boundary  rE. 
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Now,  utilizing  (11)  in  (16)  yields: 


[  'FiiidQ  +  C20  [  VFiVudQ+f2  [  'Fill dQ  +  Cg/J0  [  'FiudrE  =  C20  [  Fi^drE.  (18) 

J  Q  J  Q  J  Q  J  r £  J  Tfc 

The  weak  form  of  the  problem  is  then:  find  u  eV  and  fa  e  j  =  1, . . .  J  -  1,  such  that  Eqs.  (17)  and  (18)  are  satisfied 

VFiGVandTjeH1^). 

Having  defined  the  problem  statement  we  can  now  introduce  the  matrix-vector  form  of  the  problem.  First,  in  each  ele¬ 
ment,  we  expand  the  solution  variables  u  and  fa  using  the  same  basis  functions  used  in  the  weak  form  as  follows: 


Np  Nb 

%  =  £  Fkuk,  4>jN  =  Y) 

k= 1  k=\ 


j  =  1,2, ...  J  —  1. 


(19) 


Here,  Np  refers  to  the  number  of  points  that  Q  is  discretized  into  and  Nb  refers  to  the  number  of  points  that  rE  is  discretized 
into.  Next,  we  substitute  this  basis  function  expansion  directly  into  the  weak  form  (17),  (18)  and  simplify,  resulting  in  the 
following  matrix  form  of  the  problem: 


Mil  +  PqCIBu  +  CqLu  +/2Mu  —  C^Bfa  , 

PjMb<j>j  =  ajMb - Lbd>j_ a  -  Mb </>,_,  +  Mb<f>]+1,  j  =  l,...J-l,  '  ' ' ' 

where  the  matrices  M,  L,  B,  Mb  and  Lb  are  obtained  from  analogous  element  arrays  via  direct  stiffness  summation,  given  by 


Ne  Ne  Nt 


e=\  e=\ 


The  expressions  for  the  element  arrays  are: 

Ml  =  jf  iiij  dQe,  LI  =  jf  Vi A, d!2e, 

Bl  =  Jr  'hVjdr?,  M\  =  jr  ViVj dre,  (22) 

*■!=/  v'v'dre, 

Jre 

where  and  re  denote,  the  part  of  £2  and  r  associated  with  element  e.  Also,  ^  and  v,  are  the  locally  defined  basis  functions 
from  which  the  global  basis  functions  ( and  t,)  are  formed.  For  quadrilateral  elements  with  spectral  order  N  as  used  in  this 
analysis,  fa  will  be  discretized  into  (N+  l)2  points,  and  v,-  into  (N  +  1)  points. 

For  the  case  of  the  infinite  channel,  the  set-up  is  identical,  except  for  the  addition  of  another  set  of  auxiliary  variables  for 
use  in  the  western  boundary.  Since  there  is  no  overlap  in  the  contribution  of  the  auxiliary  variables  to  the  total  solution  (i.e.  a 
corner  condition)  there  is  no  special  handling  required  when  solving  the  dynamic  system.  The  key  difference  between  this 
work  and  that  presented  in  [23]  are  the  numerical  schemes  utilized.  This  work  uses  high-order  basis  functions  to  discretize 
the  spatial  component  of  this  problem  with  the  goal  of  improving  accuracy  over  linear  basis  functions.  In  addition,  high-or¬ 
der  time-integrators  are  used  in  conjunction  with  the  high-order  spatial  discretization  and  boundary  conditions. 


5.  Solution  of  the  dynamic  system 


The  system  described  by  (20)  and  (22)  constitute  ]  coupled  systems  of  ODEs  that  must  be  solved  for  u(x,y,t).  Our 
approach  uses  standard  kth  order  Runge-Kutta  (RK)  methods  to  integrate  the  system  in  time.  Recall  that  RK  is  used  to  solve 
first  order  ODEs,  and  as  such,  the  second  order  systems  described  must  be  converted  into  a  larger  system  of  first  order  ODEs. 
Using  the  substitution  v  =  li,  v  =  ii  and  v(j>  =  fa  is#  =  <j>  yields  the  first  order  systems: 


I 
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~ii~ 
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■  0  ■ 

0 
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Ur  _ 

i 

where  0  =  {^ ,  fa , . . . ,  fa^  }T,  a  is  the  matrix  of  assembled  terms  of  fa  A  is  the  matrix  of  assembled  terms  of  4>,  p  is  the  matrix 
of  assembled  terms  of  </>  and  ur  is  the  vector  of  terms  associated  with  fa. 

This  system  was  solved  using  a  two-stage  approach  at  each  time  step.  First,  the  auxiliary  system  was  solved  to  find  the 
component  fa  required  for  the  main  system.  Then  the  main  system  was  solved  to  find  u  at  the  next  time  step.  This  process 
was  continued  until  tf  was  reached. 
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Since  this  work  centers  on  the  use  of  high-order  spectral  elements  along  with  high-order  boundary  conditions,  high-order 
time-integrators  were  also  explored  to  examine  all  of  the  limiting  factors  of  high-order  accuracy.  We  propose  that  in  order  to 
see  the  full  improvements  of  high-order  boundary  conditions  requires  a  balance  of  truncation  errors  between  all  of  the  com¬ 
ponents  of  the  numerical  model;  this  includes  the  boundary  conditions,  the  spatial  discretization  method,  and  the  time-inte¬ 
grators.  Experiments  were  conducted  using  boundary  conditions  up  to  order  10,  SE  polynomials  up  to  order  16,  and  time- 
integrators  up  to  order  10. 

6.  Numerical  experiments 

Several  numerical  experiments  were  performed  to  solve  both  the  wave  equation  (f  =  0)  and  the  Klein-Gordon  equation. 
In  the  latter  we  used /  =  1.  Two  different  initial  conditions  were  tested.  The  first  is  a  two-dimensional  cosine  pulse  where  b  is 
the  height  of  the  channel  (in  our  case  we  used  b  =  3).  For  each  experiment,  rw  was  introduced  at  xw  =  -3.  The  initial  con¬ 
dition  is  given  by 


u(x,y,  0)  =  e-10*2  cos  •  (24) 

The  other  initial  condition  is  a  two-dimensional  Gaussian  centered  at  x  =  0,  y  =  b/ 2,  given  by 

u(x,y,0)  =  e-10x2-10(H)\  (25) 

In  order  to  see  the  effect  of  the  NRBC,  we  have  compared  our  solution  to  the  same  one  on  a  larger  domain,  i.e.  -3  <  x  <  9 
for  the  semi-infinite  channel  and  -9  <  x  <  9  for  the  infinite  channel.  Then  we  solved  the  problem  for  time  t  =  6.5  when  the 
disturbance  reached  the  wall  at  x  =  9.  We  consider  this  reference  solution  “Exact”  when  using  8th  order  basis  functions  on  a 
24  x  8  -  element  grid  (12,545  interpolation  points).  Time  integration  was  performed  with  a  4th  order  Runge-Kutta  scheme 
(RK4)  using  a  time-step  chosen  to  ensure  a  Courant  number  of  0.1,  where  the  Courant  number  is  conservatively  defined: 

Courant  number  =  C0Af -  (26) 

Y  (A*)2  +  (Ay)2 

Here,  Ax  and  Ay  are  chosen  as  the  minimum  distance  between  any  two  points  in  the  x-  or  y-directions  respectively.  This 
choice  is  made  since  the  interpolation  points  are  not  uniformly  distributed  when  using  spectral  elements.  Additionally, 
all  NRBC  parameters  C,  were  chosen  simply  as: 

Q  =  C0  =  1 . 

It  should  be  noted  that  van  Joolen  et  al.  [28]  show  that  this  choice  of  C,  is  guaranteed  to  reduce  spurious  reflection  as  the 
order  of  the  NRBC  (J)  increases. 

In  order  to  quantify  the  errors  observed  between  the  reference  and  NRBC  solutions,  we  use  the  normalized  L2Q  error  de¬ 
fined  as  follows: 


error  f2  - 


£&(«! 


,num  i,1 
k  ui 


\  ££,«) 


where  Np  is  number  of  points  in  Q  and  unum  and  uref  are  the  numerical  and  reference  solutions. 


(27) 


6.1.  Semi-infinite  channel 

For  the  first  experiment,  the  KGE  is  examined  on  a  semi-infinite  channel  with  the  NRBC  at  xE  =  3.  In  Fig.  3  we  plot  the 
exact  solution  of  the  two-dimensional  cosine  pulse  on  the  top  panel  and  the  solution  on  the  truncated  domain  with  NRBC 
on  the  bottom  panel.  The  solution  on  the  truncated  domain  with  NRBC  uses  4th  order  basis  functions  on  a  16  x  8  -  element 
grid  using  NRBC  order 7  =  4.  Qualitatively  speaking,  the  results  appear  to  show  very  little  reflection  using  the  combination  of 
high-order  basis  functions  and  high  order  NRBC  scheme. 

Representative  quantitative  results  can  be  observed  in  Fig.  4  showing  the  error  on  Q  between  the  reference  and  NRBC 
solutions  for  increasing  NRBC  order  and  various  order  basis  functions.  It  should  be  noted  that  similar  results  are  observed 
when  examining  both  the  wave  equation  and  KGE  using  either  initial  condition.  Two  main  points  can  be  extracted  from  these 
results. 

First,  it  is  clear  that  increasing  the  NRBC  order  yields  significant  gains  in  accuracy,  but  by  order  5,  the  gains  observed  by 
lower  order  basis  functions  do  not  continue.  This  may  imply  that  the  error  incurred  by  the  interior  discretization  is  more 
significant  than  the  reflection  caused  by  the  NRBC.  Second,  it  is  evident  that  by  increasing  the  order  of  the  basis  functions 
(spectral  elements),  the  gains  observed  over  simply  using  linear  (finite)  elements  are  quite  significant.  In  order  to  see  the 
gains  of  a  high-order  boundary  scheme,  one  must  have  a  high-order  interior  discretization  as  well. 
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Exact  Solution  on  Expanded  Domain 


X 


Solution  on  Truncated  Domain 


(NRBC  Order  =  4) 


-3-2-1  0  1  2  3 


X 

Fig.  3.  Dispersive  wave  guide  problem  comparing  exact  and  NRBC  solutions. 


I_Q  Error  on  Semi-Infinite  Channel  vs.  NRBC  Order 


Fig.  4.  ]}Q  error  for  dispersive  cosine  pulse  IC  using  various  order  spectral  elements. 


6.2.  Infinite  channel 

For  the  second  experiment,  the  domain  is  an  infinite  channel  with  the  NRBCs  at  x  =  ±3.  In  Fig.  5  we  plot  the  reference 
solution  of  the  two-dimensional  Gaussian  on  the  top  panel  and  the  solution  on  the  truncated  domain  with  NRBC  on  the  bot¬ 
tom  panel.  The  solution  on  the  truncated  domain  with  NRBC  uses  4th  order  basis  functions  on  a  16  x  8  -  element  grid  using 
NRBC  order  J  =  4  .  Again,  qualitatively  speaking,  the  results  appear  to  show  very  little  reflection  using  the  combination  of 
high-order  basis  functions  and  high  order  NRBC  scheme. 

Representative  quantitative  results  can  be  observed  in  Fig.  6  showing  the  error  on  Q  between  the  reference  and  NRBC 
solutions  for  increasing  NRBC  order  and  various  order  basis  functions.  This  plot  shows  results  for  the  KGE.  It  should  be  noted 
that  similar  results  are  observed  when  examining  both  the  wave  equation  and  KGE  using  either  initial  condition.  Again,  the 
conclusions  are  similar  to  those  observed  in  the  semi-infinite  channel. 
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Exact  Solution  on  Expanded  Domain 
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Fig.  5.  Dispersive  wave  guide  problem  comparing  exact  and  NRBC  solutions. 


Solution  on  Truncated  Domain 
(NRBC  Order  =  4) 
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Error  on  Infinite  Channel  vs.  NRBC  Order 


Fig.  6.  L2q  error  for  Gaussian  pulse  IC  using  various  order  spectral  elements. 


6.3.  Effects  of  time  integration  technique 

At  the  outset  of  this  work,  it  was  believed  that  at  some  point  the  improvements  realized  by  improving  the  spatial  discret¬ 
ization  and  increasing  the  order  of  the  NRBC  would  eventually  be  limited  by  the  time  integration  scheme  [29].  To  this  end, 
the  order  of  the  time  integration  scheme  was  varied  to  examine  the  effects  of  time  integration  on  accuracy  of  the  solution.  As 
has  already  been  presented,  gains  made  by  increasing  the  order  of  the  NRBC  halt  for  lower  order  spectral  elements  after 
J  =  5.  Even  for  high-order  (orders  8  and  16)  spectral  elements,  the  gains  made  by  increasing  the  order  of  the  NRBC  are  limited 
at  some  point  using  RK4. 

For  this  experiment,  consider  the  KGE  (3)  on  a  semi-infinite  domain  with  uFw  =  0.  To  ensure  that  any  boundary  or  time 
effects  are  not  masked  by  the  interior  discretization,  24th  order  spectral  elements  are  used  on  a  fine  mesh  consisting  of  4753 
global  points.  The  Gaussian  initial  condition  is  used  and  is  evaluated  until  t  =  4.  The  reference  solution  in  this  case  was  com¬ 
puted  as  described  previously,  except  this  time  using  24th  order  spectral  elements  on  a  fine  mesh  consisting  of  9457  global 
points.  Time  integration  was  performed  with  a  10th  order  Runge-Kutta  scheme  (RK10)  using  a  time-step  chosen  to  ensure  a 
Courant  number  of  0.1. 

As  can  be  observed  in  Fig.  7,  gains  made  by  improving  the  time  integration  matter  only  if  combined  with  high-order  treat¬ 
ment  of  the  boundary.  Conversely  -  gains  using  high-order  treatment  of  the  boundary  can  only  be  realized  if  there  is  a  high- 
order  treatment  of  the  time  integration. 
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Error  on  Semi-Infinite  Channel  vs.  Runge-Kutta  Order 


Fig.  7.  Error  in  the  SE  solution  of  KGE  using  NRBCs  of  various  order  as  a  function  of  time  integration  order. 

It  should  be  noted  that  these  results  (error  on  the  order  of  1CT10)  cannot  be  observed  unless  high-order  treatment  of  the 
interior  also  accompanies  the  high-order  treatment  of  the  boundary  and  time.  Several  experiments  were  conducted  which 
varied  the  order  of  the  interior,  boundary  and  time  integration.  The  clear  result  was  that  without  high-order  treatment  of  all 
components  in  concert,  convergence  to  the  reference  solution  is  stalled. 

We  believe  that  results  can  be  improved  to  machine  precision  provided  all  components  (interior,  boundary  and  time) 
have  high-order  treatment.  If  high-order  treatment  of  any  of  the  three  components  is  missing,  the  high-order  treatment 
of  the  other  components  is  essentially  wasted. 


6.4.  Long  term  results 

In  many  applications,  there  is  concern  for  the  stability  and  behavior  of  the  NRBC  scheme  in  the  long  term.  In  the  final 
infinite  channel  experiment,  we  have  run  the  KGE  with  J  =  10  using  8th  order  basis  functions  (2145  global  points)  for 
t  =  1000.  As  can  be  seen  in  Fig.  8,  by  this  time  the  waves  have  exited  the  truncated  domain  and  the  ripples  seen  in  the  figures 
are  of  order  10“5.  This  figure  shows  that,  as  expected,  in  the  absence  of  sources  in  the  domain,  the  long  term  result  is  essen¬ 
tially  zero.  This  would  imply  long  term  stability  of  the  method  when  used  in  this  context. 


7.  Conclusions 

In  this  paper,  we  have  presented  a  spectral  element  solution  to  the  Klein-Gordon  dispersive  wave  equation  on  a  truncated 
infinite  domain  in  several  settings.  Using  the  G-N  Higdon-type  NRBCs,  accuracy  was  shown  to  improve  significantly  in  the 
channel  using  a  combination  of  high-order  (spectral)  elements,  high-order  NRBCs,  and  (to  a  slightly  lesser  extent)  high-order 
time  integration  techniques.  These  results  suggest  a  balanced  approach  to  dealing  with  errors  in  solving  problems  of  this  sort 


Extended  Time  Solution 


Fig.  8.  Extended  time  (t=  1000)  KGE  solution. 
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-  namely,  to  make  improvements  in  all  major  components  of  the  solution  (interior,  boundary  and  time  discretization)  -  to 
see  improved  accuracy. 

Future  work  includes  extending  this  work  to  the  KGE  with  non-zero  advection  and  to  the  full  first  order  SWE  system.  Fur¬ 
ther  work  is  also  in  progress  to  address  quarter  and  fully  open  plane  problems  which  are  complicated  due  to  the  presence  of 
“corners”  where  two  NRBCs  intersect.  Since  the  auxiliary  variable  formulation  is  itself  a  system  of  PDEs,  well-posedness  re¬ 
quires  appropriate  boundary  conditions  -  which  are  not  obvious  -  for  each  of  the  auxiliary  variables  when  employed  in  a 
“corner-type”  configuration.  Work  is  ongoing  to  further  develop  the  auxiliary  variable  formulation  that  will  maintain  stable, 
high-order  results  in  the  quarter  and  open  domain  experiments. 
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